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ABSTRACT 

Phase retrieval refers to a classical nonconvex problem of re- 
covering a signal from its Fourier magnitude measurements. 
Inspired by the compressed sensing technique, signal spar- 
sity is exploited in recent studies of phase retrieval to reduce 
the required number of measurements, known as compressive 
phase retrieval (CPR). In this paper, l\ minimization prob- 
lems are formulated for CPR to exploit the signal sparsity 
and alternating direction algorithms are presented for prob- 
lem solving. For real-valued, nonnegative image reconstruc- 
tion, the image of interest is shown to be an optimal solution 
of the formulated l\ minimization in the noise free case. Nu- 
merical simulations demonstrate that the proposed approach 
is fast, accurate and robust to measurements noises. 

1. INTRODUCTION 

Many imaging techniques reconstruct a signal from its fre- 
quency or Fourier measurements. But in practice the phase 
information of the data in the frequency domain may not 
be available to the detecting and sensing devices, e.g., in 
X-ray crystallography [1 1. It therefore arises the problem of 
recovering a signal from its Fourier magnitude only measure- 
ments, known as phase retrieval. Other applications of phase 
retrieval include optics QJ, diffraction imaging [2], astro- 
nomical imaging [3] and magnetic resonance imaging HI, to 
name just a few. Since more than one signal can result in the 
same Fourier magnitude measurements, the phase retrieval 
problem is ill-posed. For example, a solution to the problem 
can be subject to so-called global phase ambiguity including 
a constant global phase, spatial shift and conjugate inversion. 
Though such ambiguities are acceptable in practical appli- 
cations, there may exist infinitely many solutions beyond 
these trivial associates. To resolve the problem, it is shown 
in (5J that twofold oversampling in each dimension in the fre- 
quency domain almost always specifies a unique solution (up 
to global phase) for finitely supported, real-valued and non- 
negative signals while it does not work for ID signals. Due to 
the nonlinear relationship between the magnitude measure- 



ments and the signal of interest, the algorithm design remains 
a major problem in phase retrieval. Current approach to signal 
reconstruction with the oversampling technique include the 
popular iterative projection algorithms (6) and recent convex 
optimization method |7j. 

It is of great interest to reduce the required number of 
measurements in phase retrieval since every single one of the 
measurements may have a potential time and/or power cost. 
Moreover, oversampling can be inconvenient or impossible in 
applications, e.g., Bragg sampling from periodic crystalline 
structures [8|. Compressed sensing (CS) ||9) is an emerg- 
ing technique which aims at reconstructing a high dimen- 
sional signal from its low dimensional linear measurements 
under a sparsity prior. A prosperous direction is to com- 
bine CS with phase retrieval, known as compressive phase 
retrieval (CPR) [ 8 . 10 1, which reduces the required number of 
measurements by exploiting the signal sparsity. It is shown 
in 1 10 1 that a sparse 2D image can be uniquely determined 
(up to global phase) with significantly reduced, undersampled 
Fourier magnitude measurements. Modified iterative projec- 
tion algorithms are proposed in [8} [K)}[TT| for CPR in the 
noise free case. Convex optimization method is developed in 
p2) inspired by [7] with slow computing speed. The state of 
the art is generalized approximate message passing (GAMP) 
based algorithm which is introduced in p3| and shown to 
have good performance in both recovery accuracy and com- 
puting speed. 

In this paper, nonconvex li minimization problems are 
formulated for noiseless and noisy CPR problems inspired 
by existing CS results. A desirable property is shown in the 
widely studied scenario of real-valued, nonnegative image re- 
construction that the image of interest is an optimal solution 
of the formulated l\ minimization problem in the noise free 
case. First-order iterative algorithms are proposed to solve 
the formulated problems based on the alternating direction 
method (ADM) fl4][T5) in convex optimization. Numeri- 
cal simulations are provided for simulated sparse images and 
nonnegative image to demonstrate the fast, accurate and ro- 
bust performance of the proposed approach. 



The rest of the paper is organized as follows. Section 
[2] presents our problem formulations for the noiseless and 
noisy CPR problems. Section [3]presents the ADM-based al- 
gorithms for the proposed t\ minimization problems. Section 
|4]provides the simulation results and Section|5]concludes the 
paper. 

2. COMPRESSIVE PHASE RETRIEVAL VIA l Y 
MINIMIZATION 



consider the DC (zero frequency) term of the frequency data 
which is proportional to the sum of all elements of the signal 
vector. So the DC term must be nonnegative since x° ^ 0. 
It follows from l-F^ 33 ! = |J-na; | that £\ = ^\ x°. As a 
result, 
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2.1. Problem Formulations 

We study the recovery of sparse signals from their Fourier 
magnitude measurements, known as compressive phase re- 
trieval. Inspired by the recent CS technique, we propose to 
solve the following LASSO-like optimization problem 
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where x is the signal of interest, denotes the discrete 
Fourier transform (DFT) constrained on the index set il, b de- 
notes the observed (generally noisy) Fourier magnitude mea- 
surements on and A > is a regularization parameter. The 
second term in ([T]i fits the observed magnitude data while the 
l\ norm is used to promote the signal sparsity. An equivalent 
formulation of ([TJ is 



min ||cc|| , , subject to 



\Fnx\ - b\\ 2 < e, 



(2) 



where e controls the fidelity of the reconstruction to the mea- 
sured data. It is known that ([T]l and Q are equivalent with 
appropriate choices of A and e. As A — > +00 and e — > 0, both 
([1} and (|2) reduces to 



minHcc^ , subject to |.Fna;| = b 



(3) 



in the noise free case. 



2.2. Performance Guarantee 

We note that the l\ minimization problem ([3]) has been stud- 
ied in [ 10 1 . However, there have been few theoretical results 
on its performance analysis. The following theorem deals 
with a practical and widely studied scenario where the sig- 
nal of interest is nonnegative (nonnegative means real-valued 
as well in this paper). Although far from satisfactory, this re- 
sult nevertheless provides some assurance on the reliability of 
the proposed l\ minimization approach. 

Theorem 1 If x° >z and its Fourier magnitudes b = 
\J-qX°\ is observed, where £1 contains the zero frequency, 
then x° is an optimal solution to Q. 

Proof: We need to show that the inequality > 1 1 as ° 1 1 1 
holds for any x satisfying (.Fox] = |J-h:r |. To do this, 



where the last equality follows from the positivity of x°. ■ 
Theorem[T]states that a nonnegative signal has the least l\ 
norm among all possible candidates which result in the same 
Fourier magnitudes. Though our focus is on the recovery of 
sparse signals, Theorem [T] does hold for general nonnegative 
signals and hence brings new insights to the applicable scope 
of the £1 optimization. It is noted that the proposed l\ opti- 
mization methods do not result in a unique solution. In fact, 
they suffer from at least the same ambiguities as existing ap- 
proaches to conventional phase retrieval. For example, an op- 
timal solution to either one of ([l}-([3]) after a modification of 
constant global phase, spatial shift and/or conjugate inversion 
remains optimal. 

3. ADM FOR COMPRESSIVE PHASE RETRIEVAL 

3.1. Preliminary: Alternating Direction Method 

The augmented Lagrangian alternating direction method 
(ADM) solves the structured optimization problem 

min / (x) + g (z) , subject to Ax + Bz = c, (5) 

x.z 

where / (a;) and g (z) are convex functions of x and 2, re- 
spectively. The augmented Lagrangian function of the prob- 
lem is given by 

C(x,z,y) 

= f(x)+ g (z) + (y, Ax + Bz c) 



-\\Ax + Bz-c\\ 2 2 



= f(x)+g{z) + 



a 



Ax + Bz - c 



(6) 



1 

2a 



11 2 



where y is a Lagrangian multiplier and a > is a penalty 
parameter. Starting with y° and z°, the ADM iterates as fol- 
lows: 



x k+1 — argmin£ (x,z k ,y k ) 
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where /3 € ^0, v/ ^, +1 ^ guarantees the convergence under 

some technical assumptions fl4| . The ADM is very efficient 
when explicit solutions are available for |7]i and (|SJ. The 
ADM has been a popular approach to solutions of large scale 
problems, since it can typically produce a modestly accurate 
solution within a few tens of iterations though its convergence 
to high accuracy may be slow 1 1 5 1 . 

3.2. ADM for Compressive Phase Retrieval 

Problems ([T]i-(|3]l are nonconvex due to the presence of the 
magnitude operator. In this section we solve ([T])-(|3]l using the 
ADM. Problem ([T| can be formulated into 



^\\\zn\-b\\l 



, subject to Fx — z = 0. 



(10) 



According to (|5j, we see in this case that / (x) = 
convex while g (z) — ^ \\\zq\ — b\\ 2 is nonconvex. By rec- 
ognizing that the DFT is unitary, x can be explicitly updated 



using a soft thresholding operator (later given in ( 12 1). Let 
s = Fx + —y. Then the augmented Lagrangian involving z 
is expressed as 
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where the '=' in the inequality holds if sgn (zi) = sgn (.s^) 
for i £ SI and Zj = Si for i f2, and C is a constant indepen- 
dent of z. As a result, the augmented Lagrangian obtaines the 
minimum at the z such that Zi — Ab ^+^l s 'l sgn (sj) for i e Q 
and Zi = Si for i ^ SI. So the alternating direction algorithm 
for ([T| is summarized as follows: 
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otherwise, 
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(14) 



where s fc+1 = Fx k+1 + ^y k , S a (w) — sgn (w) ■ (\w\ — A) + 
is a soft thresholding operator with (■)+ = max(-, 0). 
Next, Q can be formulated as 

min \\x\l-. + Is {z) , subject to Fx — z = 0, (15) 

where S = {z : \\\zq\ — b\\ 2 < e}, and Ig is an indicator 
function with Ig (z) = if z £ S, and Ig (z) — +oo other- 
wise. As a result, we have / (x) = , and g (z) = Ig (z) 



is nonconvex since the set S is nonconvex. Denote Vg the 
projection (in Euclidean norm) onto S. According to the 
ADM we obtain the algorithm for Q which is the same as that 



for |lji except for the update rule of z, z 
So, provided 6^ Owe have 

yk+ i_j(e\ s k+1 \ + (i-9)b t ) S gms 
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otherwise, 
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1 is defined as 



before. A detailed derivation of ( fTo*) is omitted. 

AsA^ +oo and e -> 0, both ( pj) and ( |T6] > reduce to the 
following update rule of z for solving ([3]): 



,fc+i 



if i G O, 
otherwise, 



(17) 



while again the updates for x and y are the same as (fT~2|> and 



14 1, respectively. 

Unfortunately, the proposed ADM algorithms cannot be 
guaranteed to converge. In fact, the convergence issue re- 
mains a difficult problem in the nonconvex phase retrieval 
problem and a mathematically rigorous analysis has not been 
found for many existing algorithms for conventional phase re- 
trieval including HIO, HPR and RAAR [6]. The performance 
of the proposed algorithms will be numerically studied in Sec- 
tion]?] where we observe that the proposed algorithms can 
either converge to a good result or not converge at all. 

4. NUMERICAL SIMULATIONS 

We provide numerical simulations to demonstrate the per- 
formance of the proposed ADM algorithms in this section. 
2D Random sparse images and a nonnegative satellite image 
shown in Fig. 3(a) are considered. Due to the page limit, a 



thorough comparison of our method with existing approaches 
will be reported in a future publication. 

Algorithm implementation: The algorithms start with 
y° = and a random z° such that \zq\ = b and Zj = for 
i ^ Q. a is set such that a fixed portion, p, of the elements 
of x 1 are kept nonzero with p = 0.6. We set f3 = 0.5 for 
random sparse images and (3 — 0.8 for the satellite image. 

max(a||z fc -z fc - 1 || .Fx k -z k ) 

The algorithm is terminated if — ^ 2 < 

10 -3 or a maximum number of iterations, 500, is reached. 

In the first simulation, we consider random images in the 
noiseless case and study the success rate of recovery versus 
the sparsity level K. Three types of images are considered in- 
cluding complex, real and nonnegative ones to test whether 
they have different recovery performance since Theorem [T] 
holds only for nonnegative signals. We consider images of 
dimension 16 x 16 (number of image pixels N = 16 2 = 256) 
with image sparsity level K varying from 1 to 35 and acquire 
a number of M = A^/2 = 128 random Fourier magnitude 
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Fig. 1. Success rates of recovering complex, real and nonneg- 
ative random images with respect to the sparsity level K. The 
number of image pixels and Fourier magnitude measurements 
are N = 256 and M — 128, respectively. 



measurements. An image is claimed to be correctly recovered 

if the relative root mean squared error (RMSE) ^ f N ^ 2 < 

\\ x 1I2 

10~ 2 , where x denotes the recovered image after removing 
the possible ambiguities of constant global phase, spatial shift 
and conjugate inversion. For each combination of K and the 
image type, 200 trials are repeated. The success rates of re- 
covery are presented in Fig. [T] In the case of a moderate 
or low sparsity level, the proposed algorithm is observed to 
converge to the correct solution in most trials starting from a 
random point. Moreover, it is shown that nonnegative images 
have the highest success rate of recovery followed by real and 
then complex ones. The inherent reason will be explored in 
future studies. 

The second simulation studies the variation of the recov- 
ery error with respect to the noise energy. We fix the sparsity 
level K = 8 and set N and M as before. Complex random 
images are scaled to unit Frobenius norm. After obtaining the 
noiseless Fourier magnitude measurements, a Gaussian ran- 
dom noise is added with the noise energy e = || |.Ffia; | — b\\ 2 
varying from to 0.2 with a step size of 0.005. The aver- 
aged relative RMSE is presented in Fig. |2j where it is shown 
that the recovery error grows approximately linearly with the 
noise level e. 

The last simulation studies the reconstruction of a nonneg- 
ative 256 x 256 satellite image (N = 256 2 = 65536) which 



is shown in Fig. 3(a) The satellite image has a sparsity ratio 
K/N ps 0.14. The number of Fourier magnitude measure- 
ments, M, is set such that M/N = 0.2, 0, 5 and 1, respec- 
tively. In each setting, a Gaussian random noise is added to 
the noiseless magnitude measurements such that the signal 
to noise ratio (SNR) is 30 dB. A fixed number of iterations, 
200, is used to reconstruct the images. Our simulation results 



Fig. 2. Reconstruction error of random sparse images vs. the 
noise level e with N = 256, M = 128 and K = 8. 



are presented in Fig. [3] where, remarkably, a faithful recon- 
struction is obtained with only 20% of the Fourier magnitude 
measurements. Without accounting for the effects of the pos- 
sible ambiguities mentioned before, the three reconstructed 
images have SNRs of 11.6, 15.3 and 18.2 dB, respectively. 
Moreover, the proposed algorithm is very fast and takes 7.8, 
8.8 and 9.3 s, respectively, to obtain the reconstructed images 
using Matlab v7.7.0 on a PC. 

5. CONCLUSION 

The noiseless and noisy compressive phase retrieval problem 
was studied in this paper while the l\ norm is exploited to 
promote the signal sparsity inspired by compressive sensing. 
The optimality of the formulated l\ minimization problem 
was proven for nonnegative signals in the sense that the signal 
of interest is an optimal solution of the i\ minimization prob- 
lem in the noise free case. Efficient alternating direction algo- 
rithms were proposed for the problem solving and promising 
results were presented to demonstrate their performance. 
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